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1.  INTRODUCTION 


~  -- '  Equations  for  the  motion  of  a  continuous  medium  capable  of  supporting 

shear  stresses  are  reported  by  a  number  of  authors,  but  the  complete 

set  does  not  appear  to  have  been  published  in  a  form  suitable  for  solution 
by  Eulerian  hydrodynamic  codes.  Specifically,  the  Eulerian  equations  of 
motion  for  a  compressible  medium  acted  upon  by  a  general  stress  tensor  are 
required.  In  this  volume  the  equations  of  motion  sure  discussed  starting 
from  the  principles  of  mass,  momentum  and  energy  conservation  for  finite 
masses.  From  these  the  differential  equations  of  motion  are  derived.  The 
constitutive  equations  relating  stress  and  strain  are  required  to  complete 
the  mathematical  description. 

It  is  also  shown  that  the  difference  equations  for  hydrodynamic  codes 
can  be  conveniently  obtained  from  the  integral  form  of  the  conservation  laws. 
This  is  more  convenient  for  a  parti cular\ coordinate  system  than  using  the 
general  method  of  tensor  analysis  and  covariant  differentiation,  which  tends 
to  be  awkward  when  deriving  equations  of  motion  in  a  specific  curvilinear 
coordinate  system. 

The  nature  of  a  medium  is  specified  by  its  equation  of  state,  which 
is  vised  to  calculate  the  pressure  from  the  density  and  specific  internal 
energy,  and  u  tensor  constitutive  equation  relating  deviator  stresses,  de¬ 
viator  strains  and  their  rates.  Constitutive  equations  for  the  shear  stresses 
are  discussed  and  the  appropriate  forms  for  elastic,  elastic-plastic  and 
rigid-plastic  solids  as  veil  as  for  viscous  fluids  are  presented.  ~v  — 
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2.  A  USEFUL  RELATION 

Before  proceeding  to  discuss  the  specific  conservation  lavs,  it  will 
prove  useful  to  formulate  an  integral  relation  connecting  the  Eulerian  and 
Lagrangian  descriptions  of  fluid  motion.  This  relation 


a_ 

dt 


Eds 


(1) 


expresses  the  physical  notion  that  the  rate  of  change  of  an  integral  over 
a  volume  fixed  in  space  (7)  is  the  difference  of  two  terms,  a  Lagrangian 
derivative  (taken  with  respect  to  the  moving  element  of  mass  contained  in 
V)  and  a  rate  associated  with  mass  transport  out  of  V. 

To  prove  this  relation  consider  the  volume  integral 


J(t) 


'/v 


H(t)  dv 


over  the  volume  V  containing  a  fixed  element  of  mass.  In  an  Interval  At 
the  change  of  J  can  be  written 


A  J 


*  fi  H(t+At)  dv  -  f  H(t)  dv 

Jv(  t+At)  yv(t) 


which,  to  the  first  order  in  At,  can  be  expressed  as 


J  *  f  _  H(t)  dv  +  At  f 

J  v(t+At)-v(t)  yv(t) 


H(t)  dv 


By  reference  to  Fig.  \  (where  S  is  the  surface  bounding  V,  n^  denotes 

the  i  component  of  the  outward  drawn  normal,  and  the  u^  are  the  components 

of  the  velocity  vector)  it  can  be  observed  that  the  thickness  of  the  volume 
—  * 
swept  out  by  S  in  time  At  is  u,n^  At. 


Tensor  notation  is  used  throughout  this  report,  so  that  in  this  ex¬ 
pression  summation  over  repeated  indices  is  to  be  understood.  The  text  by 
Brillouin'^'  gives  a  good  discussion  of  the  physics  and  mathematics  appro¬ 
priate  here. 
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V(t+At  i 
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Fig.  1 — Illustrating  the  motion  of  the  volume 
V  and  the  volume  it  sweeps  out 


Then  the  expression  above  for  AJ  becomes 


AJ  •  f  H(t)(u  n.At)  dS 


+  tt/v 


H(t)  dv 


If.  in  the  integral  over  V  shown  above,  V  is  replaced  by  its  instantaneous 
counterpart,  V,  which  is  fixed  in  space,  then  the  time  derivative  can  be 
taken  outside  the  integral.  When  that  is  dooe  and  the  resulting  equation 
is  divided  by  At,  the  expression  for  the  change  of  J  which  was  written 
as  Eq.  1  is  obtained: 

ft fv  H(t)  4v  ■  it  £ H(t)  dT  -f-s  Vi Bds  • 


denotes  the  derivative  of  J  as  the  mass  in  V  moves  along  with  the  fluid. 
Special  cases  of  interest  are  those  in  which  the  integral  on  the  left  of 
Eq.  1  denotes  the  mass,  momentum  and  energy  within  a  fixed  element  of  volume, 
and  will  be  discussed  in  the  next  paragraphs. 


3.  CONSERVATION  OF  MASS 

If  we  put  for  H  the  density,  p,  of  the  fluid,  then  the  first  in¬ 
tegral  on  the  right  of  Eq.  1  bee  ones  the  change  of  mass  in  V,  which  vanishes 
since  V  is  defined  to  contain  a  fixed  element  cf  mass.  Thus  we  find 


» dv  ■  -fa c  v» 


dv 


(2) 


and  applying  Gauss'  theorem 

fs  Fi  ni ds  m/v  Fi,i  * 

we  find 

P  +  (Pui)#i3  dv  »  0 

_  \ 

and,  since  the  volume  V  is  arbitrary, 

P  +  (Pui)^i  *  0 

A  compact  notation  is  obtained  with  the  convective  derivative: 


-  a-  +yv  -a- 

Dt  at  Za  i  axA 


The  final  form  of  the  mass  equation  is,  t*w  current  notation: 


(3) 


(*) 


(5) 


2>  + 


Dt  PUi,i  “  ° 


(6) 


i'.  CONSERVATION  OF  MOMENTIM 

To  calculate  the  momentum  of  a  finite  volume,  we  require  unit  vectors 
having  a  fixed  direction,  D,  in  space.  The  components  of  these,  vectors  may, 
of  course,  vary  from  point  to  point  in  a  curvilinear  coordinate  system. 


4 


Letting  b^  denote  the  components  of  such  a  vector,  the  D- component  of 
momentum  of  the  mass,  in  V  is  given  by 

Jv  "ui  ”i dv 

The  rate  of  change  of  momentum,  by  Newton's  law,  is  equal  to  the  component 
in  the  fixed  direction,  D,  of  the  surface  tractions  over  If,  i.e. 


St  Jv  "ui  bt  dv  -fg  bi  “j 


dS 


where  is  the  total  stress  tensor,  the  force  per  unit  aree  in  7>  is 

n^  *  F^  and  the  component  of  force  in  the  direction  of  b^  is  F^  b^ 


n^  b^  The  general  relation  of  Eq.  1  then  becomes 


&  i  bi dv '  /  Tu b 

✓  v  c 


V  ~  s 

where  the  two  integrals  have  been  combined  by  putting 

Tij  *  °u  *  "Vi  uj 

Applying  Gauss '  theorem  we  find,  as  before, 

a 


St  (Pui  V  *  (TiJ  \\i 


(7) 


(8) 


(9) 


which  is  similar  to  Eq.  X.96  of  Ref.  2.  In  rectangular  coordinates  this 
becomes,  after  combining  with  Eq.  6 

Du. 

■  (10) 


5.  CONSERVATION  OF  ENERGY 

The  energy  equation  is  obtained  by  a  process  similar  to  those 
followed  above.  The  total  energy  per  unit  mass  is 

E  -  I  +  fcuj  Uj  (11) 
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where  I  denotes  the  internal  energy  per  unit  mass-.  Since  the  rate  of 
change  of  the  total  energy  within  V  must  equal  the  rate  at  which  work  is 
done  by  the  surface  tractions  we  can  write 

Ft »(1  *  *uj  uj>  dv  -f-s  °i}  UJ  ni  48  (12) 

Putting  pE  for  H  in  Eq.  1  we  obtain  the  energy  equation 

ft fv  pE  dT  mf%  ui  “j 48  (13) 


which,  by  the  same  arguments  as  before,  leads  to  the  result 

ft  (',E)  ’  3^  [cijuj-PEui5  (l4) 

This  can  be  combined  with  the  momentum  equation  to  obtain 


'  DI 

p  Dt  "  aij  Uj,i 

which  is  just  the  first  law  of  thermodynamics  in  hydrodynamic  notation. 

\ 

6.  THE  FINITE  DIFFERENCE  EQUATIONS 


(15) 


The  integral  relations  of  Eqs.  2,  7  and  13  can  be  used  to  write  out 
directly  the  flow  equations  in  finite  difference  form.  This  approach  seems 
natural  since  it  suggests  a  finite  difference  scheme  in  which  the  conservation 
f  mass,  momentum  and  energy  are  naturally  assured,  and  a  physical  interpre¬ 
tation  of  each  step  in  the  iJumerical  computation  is  possible.  An  appropriate 
control  voluca  for  flow  with  axial  symmetry  is  sketched  in  Fig.  2.  It 
consists  of  half  a  toroid  of  rectangular  section.  The  choice  of  a  half  ring 
allows  for  a  rigorous  determination  of  the  momentum  equation  for  radial  velo¬ 
cities,  but  for  the  other  equations  a  complete  toroid  would  do  just  as  well. 

If  the  integrands  in  Eq.  2  for  mass  conservation  are  taken  to  be 
constant  and  V  is  the  half-ring  of  Fig.  2,  then  the  surface  integral  vani¬ 
shes  over  S3  an.d  in  view  of  the  zero  angular  velocity,  but  the  in¬ 
tegral  over  the  other  four  surfaces  is  finite. 
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Fig.  2 — The  shape  of  an  element  of  volume,  V,  appiopriate  for 
deriving  the  Euleilan  finite  difference  equations  in  cylindrl 
cal  coordinates  is  illustrated  in  the  sketch. 

Then  Eq.  2  becomes 


...  2  6  > 

p  V  «  -  2«  Jpru  |  Az  +  pv  |  (r+^Ar)  Ar>  (16 

l  1  5  > 

In  this  and  subsequent  equations  we  specialize  the  notation,  u^,  in  view  of 
the  cylindrical  geometry,  to 

-  u,  u2  -  0,  u^  ■  v  (D 

Eq.  16  is  the  mass  transport  equation  used  in  the  finite  difference  scheme, 
except  for  replacement  of  the  time  derivative  by  a  .finite  difference. 

The  radial  momentum  equation  is  obtained  by  means  of  Eq.  7  in  which 
b^  is  the  unit  vector  in  the  (fixed)  x  direction.  Then,  with  0  the 
azimuth  angle  to  the  x-axis, 

b^  *  cos  0,  bg  ■  -  sin0,  b^  ■  0  (l l 

In  view  of  the  symmetry  several  components  of  T^  vanish,  viz. 


and  the  integral  equation  is 

rr  C  p  u  cos  0  d v  *  /*  (T.  .  cos  0  -  T0.  ^in  0)  n.  dS 

dVv  J  s  1J  J 

As  before,  it  is  assumed  that  p  and  u  are  constant  in  V  and  T.^  is 

constant  on  each  of  the  faces  of  the  cell.  Carrying  out  the  integrals,  we 

find,  upon  multiplication  by  n, 

a  /  2  6  . 

(pu?)  »  2 n  <  rT^  |  Az  -  T^  Ar  Az  +  T^  {  Ar(r+$Ar)j  (20  ) 

i  *  5 

For  the  momentum  equation  in  the  axial  direction  we  put  for  h^ 
b..  »  0,  b9  ■  0,  b-  ■  1 


8 


which,  with  Eq.  7>  leads  to 


.  ,26. 

(pv  V)  »  2k  |rT31  Lz  +  |  (r+£dr)  h  r j-  (21 

The  finite-difference  energy  equation  follows  by  carrying  out  the 
integrals  of  Eq.  13  over  V,  leading  to 

a  i  2 

—  (E yf)  «  2k  jrto^  -  pE)  u  +  cr31  v]  |  A z 

6  > 

+  [o13  U  +  (g33  -  pE)  v  j  (r+|Ar)  ArJ.  .  (22 

In  sunroary,  we  have  four  equatione  for  the  conservation  of  mass,  two 
components  of  momentum  and  the  energy.  By  means  of  Eq.  27  we  define 
through  A^,  and  in  addition  we  put  m  for  pV  to  obtain  the  equations  in 
their  most  compact  form.  These  equations  are  essentially  those  vised  in  the 
computer  program,  except  for  the  finite  time  step. 

|f-  -  O^-pvAK  (2; 


d  (mu) 


if- ^ 


a(mv) 


ir-V^V, 


2&L  C(ou  -  PE)  u  ♦  v]  *1* 


+  [o,,  U  +  (a,,  -  pE)  v]  A| 


The  surfaces  referred  to  are  illustrated  in  Fig.  2,  and  the  value  of  A  on 
the  i—  surface  ia  denoted  by  A^. 


K 


A1  «*  2«  r  Az,  Ag  »  2n  (r+Ar)  Az 

A3  “  \  m  2n  Ar  Az 
A^  *  A^  ■  2n  (r+^Ar)  Ar 


1 


(27) 


7.  DIFFERENTIAL  EQUATIONS  FOR  FLCW  WITH  CYLINDRICAL  SYMMETRY 

For  completeness  the  flow  equations  in  differential  form  which  are 
obtained  by  passing  to  the  limit  of  vanishing  Ar  and  Az  are  noted,  below. 


ap  ,  _  i  MehI  .  M 
at  r  ar  az 


(26) 


1  a(Tnr)  i  8T«  Tgg 

at  r  •  ar  az  r 


a(pv)  -  1  a(  t  aT33 
at  r  ar  az 


(29) 

(30) 


alfiE] 

at 


P^{c(aii"  pE)  u  +  °3i  r} 

*  {a13  u  *  (a33"pE)  V} 


(31) 


8.  ISOTROPIC  AND  DEVIATOR  COMPONENTS  OF  STRESS  AND  STRAIN 

The  stresses  a.  ,  can  be  separated  into  two  parts,  a  pressure  (iso- 
tropic)  component  p,  which  is  the  average  of  the  principal  stresses,  and  a 
deviator  component,  S^,  which  results  from  subtracting  out  the  pressure 
term  from  the  total  stress.  Thus  we  can  write 


°ij  '  •  p  h)  *  su 

(32) 

and  it  follows  that 

*  -  -  Ki 

(33) 

and 

. 

S  »  0 
°ii 

(3k) 
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The  negative  sign  preceeding  the  pressure  in  Eq.  32  is  required  to  conform 
with  the  visual  conventions  that  pressure  is  positive  for  material  in  com¬ 
pression  and  stresses  are  positive  for  material  in  tension. 

In  a  similar  manner  the  velocity  gradient  can  he  resolved  into  a 
rotation,  a  compression  and  a  distortion.  To  accomplish  this  the  stmin- 
rate  tensor  c .  ^  is  defined  as  the  symmetric  component  of  the  gradient  of 
the  velocity  field 


cu 

and  the  skew  part 

“u 


uJ,i) 

(35) 

(36) 

of  the  tensor  u.  .  defines  the  vorticity.  The  commas  in  these  expressions 
i*  J 

denote  differentiation  in  the  sense  of  tensor  analysis.  Since  the  flows  of 
interest  in  this  report  have .axial  symmetry,  the  matrix  of  strain- rates  is 
noted  helcw  for  that  special. case: 


Bu 

dr 


*<£♦£) 


<«u>  ’  < 


u 

r 


+  £>  0 


Bz 


(37) 


The  strain  rate  can,  in  turn,  he  broken  up  into  an  isotropic  part  involving 
the  rate  of  change  of  density  and  a  deviator  part  obtained  hy  subtracting 
out  this  isotropic  ''omponent.  Specifically, 


'ij 


ciJ  "  ^*kk  6iJ 


(38) 


A  relation  between  strain-rate  and  density  is  obtained  by  combining  Eqs.  6 
and  35: 


1  Do  1  DV 
8ii-  Ui,i  “•  p  Dt  VDt 


(39) 
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The  pressure  is  determined  by  the  equation  of  state,  p(E, p).  A  general 

discussion  of  current  information  on  the  equation  of  state  is  given  by 

Brush  and  a  specific  empirical  fit  to  the  existing  data  is  described 
(4) 

by  Tillotson.  ' 

To  complete  the  description  of  the  flew  a  specific  relation  between 
stress  and  straija  ^s  needed.  A  general  form  of  such  an  equation  is  displayed 
and  discussed  below  to  place  the  particular  equations  used  to  formulate 
the  OIL-RPM  code  in  context.  Subsequently,  it  will  be  specialized  to  des¬ 
cribe  the  rigid-plastic  model  used  in  the  computer  program. 


An  equation  of  the  form 


JL  s 

^  1,5 


+  b  S 


1.5 


(40) 


relating  strain-rate,  stress-rate  and  stress  includes  the  description  of 
classical  viscosity,  elasticity,  end  the  Prandtl-Reuss  equations  of  plastic 
flow  as  special  cases.  Here  b  and  p,  are  scalars,  but  not  necessarily 
constants.  The  dots  denote  differentiation  following  a  set  of  axes  fixed 
in  the  material.  Thus,  special  provision. has  to  be  made  in  the  derivatives 
if  there  is  rotation  or  translation  of  the  material. 

If  p,  *  •  and  b  «  l/2pv  the  equations  describe  classical  viscosity 
with  v  the  kinematic  viscosity  and  pv  the  dynamic  viscosity  of  the 
material.  The  full  equations  are  then  the  Navier-Stokes  equations  if  the 
additional  requirement  is  made  that  the  density  be  held  constant,  as  in 
incompressible  flow. 


Elasticity  is  described  by  putting  b  *  0,  and  in  that  case  p,  is 
the  rigidity  modulus  of  the  material,  which  is  one  of  the  two  Lam£  constants. 

A  second  elastic  constant  is  included  in  the  equation  of  state  of  the  material, 
which  for  moderate  internal  energy,  I,  and  near-normal  density  is  approximated 
by  p  *  A(p/pQ-l),  where  A  is  the  bulk  modulus  in  the  notation  of  this 
report.  The  other  Lam£  coefficient  is  expressible  in  terms  of  A  and  p, 
by  \  «  A  -  §p. 


12 


mgmsmsms*-- 


The  Prandtl-Reuss  equations  for  an  elastic- plastic  medium  with  a 
von  Mises  yield  criterion  are  discussed  by  Hill, ^  Thomas and  numerous 
other  authors  on  plasticity.  If,  in  addition  to  Eq.  40,  we  require  the 
von  Mises  yield  criterion 

su  su  - 

where  Y  is  the  yield  stress  in  pure  shear,  we  obtain  an  explicit  expres¬ 
sion  for  b. 


which  completes  the  specification  of  the  medium,  except  for  the  possibility 
that  the  deformation  is  elastic.  In  that  case  b  ■  0.  At  each  step  in 
the  calculation  the  deformation  has  to  be  tested  to  see  whether  the  new 


stress  is  elastic,  l.e.  if 


sij  su  < 


or  plastic,  which  is  the  case  when  the  inequality  is  violated,  as  in  Eq.  4l. 

The  full  Prandtl-Reuss  equations  have  to  account  for  material  rotation 
by  means  of  additional  terms  which  were  originally  described  by  Jaumann 
(see  Prager^  for  a  discussion^  and  are  described  by  Thomas^  on  p.  88. 

In  our  notation 

su  ■  8t  su  *  skj  “ki  *  sik  “Vj  >  (M0 

where  oj. ,  is  given  by  Eq.  36  and  3.  denotes  an  ordinary  time  derivative 
in  fixed  axes. 

V 

The  rigid- plastic  model  constitutes  an  Important  simplification  which 
is  useful  when  the  strains  are  large  compared  to  the  value  at  yield,  as  for 
the  "near  field"  in  hypervelocity  impact.  In  this  approximation  the  elastic 

a 

portion  S.  ./2\i  of  the  total  strain  is  set  to  zero,  resulting  in  the  con- 
stitutive  relation 
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This  may  be  thought  of  as  an  approximation  in  which  ^  is  very  large.  In 
this  limit  of  an  infinite  modulus  of  rigidity  it  is  said  that  the  material 
is  described  by  a  "rigid- plastic"  model.  Setting  ^  equal  to  infinity 
eleminates  the  need  for  calculating  the  Jaumann  terms,  in  Eq.  44,  and  allows 
the  stresses  to  be  calculated  directly  from  the  strain-rate.  To  do  this 
Eq.  45  is  multiplied  by  itself,  with  the  result 


b2  S 


id 


Accounting  for  the  von  Mises  condition,  Eq.  4l,  we  have 


n  Y 


(46) 


\ 
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